Predictive potentials of glycosylation-related genes in glioma prognosis and their correlation with immune infiltration

Glycosylation is currently considered to be an important hallmark of cancer. However, the characterization of glycosylation-related gene sets has not been comprehensively analyzed in glioma, and the relationship between glycosylation-related genes and glioma prognosis has not been elucidated. Here, we firstly found that the glycosylation-related differentially expressed genes in glioma patients were engaged in biological functions related to glioma progression revealed by enrichment analysis. Then seven glycosylation genes (BGN, C1GALT1C1L, GALNT13, SDC1, SERPINA1, SPTBN5 and TUBA1C) associated with glioma prognosis were screened out by consensus clustering, principal component analysis, Lasso regression, and univariate and multivariate Cox regression analysis using the TCGA-GTEx database. A glycosylation-related prognostic signature was developed and validated using CGGA database data with significantly accurate prediction on glioma prognosis, which showed better capacity to predict the prognosis of glioma patients than clinicopathological factors do. GSEA enrichment analysis based on the risk score further revealed that patients in the high-risk group were involved in immune-related pathways such as cytokine signaling, inflammatory responses, and immune regulation, as well as glycan synthesis and metabolic function. Immuno-correlation analysis revealed that a variety of immune cell infiltrations, such as Macrophage, activated dendritic cell, Regulatory T cell (Treg), and Natural killer cell, were increased in the high-risk group. Moreover, functional experiments were performed to evaluate the roles of risk genes in the cell viability and cell number of glioma U87 and U251 cells, which demonstrated that silencing BGN, SDC1, SERPINA1, TUBA1C, C1GALT1C1L and SPTBN5 could inhibit the growth and viability of glioma cells. These findings strengthened the prognostic potentials of our predictive signature in glioma. In conclusion, this prognostic model composed of 7 glycosylation-related genes distinguishes well the high-risk glioma patients, which might potentially serve as caner biomarkers for disease diagnosis and treatment.


Acquisition and collation of gene expression matrix and clinicopathological data
First, gene expression sequencing of 689 glioma samples from TCGA (LGG: n = 523, GBM: n = 166) and 5 paraneoplastic samples from TCGA and 1152 normal brain tissues from GTEx database were downloaded from University of California Santa Cruz (UCSC) Data.To produce gene expression matrices for 1846 samples, the data was log2 (TPM + 1) processed and gene ID transformed.We extracted clinical case factors of glioma patients from the TCGA and cBioPortal databases (https:// www.cbiop ortal.org/), respectively, and extracted survival outcome, survival time, age, gender, tumor grade, histological subtype, 1p/19q co-deletion status, radiotherapy status, MGMT promoter methylation, and IDH mutation status.The data information of glioma patients in the Chinese Glioma Genome Atlas (CGGA) database was used as external validation, and the data information of 1018 glioma cases with dataset IDs mRNAseq_693 and mRNAseq_325 and dataset ID mRNA sequencing (nonglioma as control) of 20 non-tumor brain tissues to obtain the gene expression profile matrix.To eliminate batch effects, the data was processed using the "limma" and "sva" packages.

Differential analysis and clustering analysis of glycosylation-related regulators
Following that, we looked at the differential expression of glycosylation-related regulatory genes in glioma and normal brain tissues, using the following criteria: adj-P < 0.05 and differential multiplicity |Log2FC > 1| were regarded differently expressed.We used violin plots and expression heat maps to plot correlation analyses for genes with differential multiplicity > 2. Consensus clustering is a data mining technique for identifying a group of unknown genomes that share biological characteristics.We split 689 glioma patients into two different subgroups and used PCA master analysis using the "limma," "ConsensusClusterPlus," and "scatterplot3d" packages to confirm the reliability of the consensus clustering results.

Enrichment analysis and survival analysis
The functional annotation of differentially expressed genes between the two groupings was done using Gene Ontology (GO) and the Kyoto Encyclopedia of Genes and Genomes (KEGG) 13 , and Kaplan-Meier survival analysis was performed to assess differences in survival rates and plot survival curves.The packages "clusterProfiler" and "survival" are used to carry out the analysis phases.The variations in clinicopathological factors between the two clusters were analyzed using the chi-square test, which was shown as a heat map.

Immunocorrelation analysis and clinical correlation analysis
To investigate the variations in immune responses between the two risk groups, we used the ssGSEA algorithm to evaluate the infiltration scores of 28 immune cells and immune-related pathways in the two groups based on the "GSVA" package's risk ratings.Because immune checkpoint expression levels may be related to immunotherapy response, the sensitivity of immune checkpoints was compared between high-and low-risk groups using the "ggpubr" package.A swarm map of associations between predicted genes, risk scores, and clinicopathological factors was created using the "beeswarm" package.

Cell culturing
The U251 and U87 glioma tumor cell line, purchased from the Animal Zoology Department of Kunming Medical University, were cultured and used to generate derivatives.U251 and U87 cells were routinely grown in DMEM/high glucose (Hyclone, USA) containing 10% Foetal Bovine Serum (FBS; Hyclone, USA) and 1% penicillin-streptomycin solution (PSS, Hyclone, USA).After being washed with phosphate-buffered saline (PBS; Hyclone, USA) one or two times, the cells were digested for 2-3 min (min) with 0.25% trypsin (1-2 mL; Gibco).The cell digestion was terminated using complete medium.After centrifugation (800-1000 rpm for 5-8 min) and resuspension, the cell suspension was collected again and the cells were placed in a 25 t (3 mL) culture flask at a density of 4 × 10 5 cells/mL and placed in the incubator.

Small interfering RNA (siRNA) transfection assay
SiRNA is the most commonly used tool for gene silencing experiments, which can direct the RISC system to specifically shear the target RNA inside the cell thereby inhibiting the expression of the target gene.The siRNAs provided by RIBO are modified with 5'Chol + 2'OMe.Before transfection, riboFECT™ CP Buffer (1X) was prepared by diluting riboFECT™ CP Buffer (10X) using PBS.After removing the riboFECT™ CP Reagent, it was fully shaken in a vortex shaker and left at room temperature to return to room temperature before use.First, siRNA solutions were diluted in 1X riboFECT™ CP buffer to obtain final concentrations of 50, 75, and 100 nm.Next, RiboFECT ™ CP reagent was added to prepare the transfection complex, which was then added to the cell cultures.RNA was extracted 48 h after transfection for validation via RT-qPCR and the most efficient concentration was selected for formal experiments.The cell culture plates were incubated in a CO 2 incubator at 37 °C for 48 h (h) for siRNA silencing effect assay and cell viability assay.

Cell Counting Kit 8 (CCK8) assay
Cell viability after transfection was assayed by CCK-8.The cells were collected in the logarithmic growth phase, and the concentration of cell suspension was adjusted by adding 100 μL to each well of a 96-well plate to achieve a cell density of 3000 cells/well for the cells to be tested.The cells were incubated in 5% CO 2 , 37 °C incubator for 24 h.After transfecting the cells with different concentrations of interfering RNA, the plates were returned to the incubator and incubated for 48 h.After 10 μL CCK-8 solution was added to each well, incubation was continued for 2 h in the cell culture incubator, and the absorbance of each well was measured at 450 nm by an enzyme meter.Cell viability = [(experimental wells − blank wells)/(control wells − blank wells)] × 100%.

Statistical analysis
All statistical analyses were performed using R software (version 4.2.1) and SPSS software (version 13.0; SPSS, Inc., Chicago, IL, USA).Data between the two groups were analysed using the Student's T test.Data are presented as the mean ± standard deviation (SD).Statistical significance was set at p < 0.05.

Identification of glycosylation-related regulators expressed in patients with glioma
First, 552 glycosylation-related genes were chosen for follow-up studies by downloading glycosylation-related regulators in five gene sets from the MSigDB database.After deleting duplicate genes, a total of 552 glycosylationrelated genes were selected for follow-up investigations (Supplementary Table 1).We used the TCGA-GTEx database to compare the expression of 552 glycosylation-related genes in gliomas (n = 689) and non-tumor brain tissues (n = 1157), finding that 43.61% (n = 215) were up-regulated genes and 8.52% (n = 42) were down-regulated genes, with the top ten up-and down-regulated genes marked in the volcano plot (Fig. 1A,B).Following that, we chose 53 glycosylation-related genes in gliomas with |Log 2 FC > 2| for further investigation (Supplementary Table 2).The Pearson correlation analysis was used to look at the relationships between genes, and it found a strong positive correlation (0.94) between SPTB and SPTBN4 and a strong negative correlation (− 0.53) between ST8SIA3 and MAGT1 (Fig. 1C).We can see that there are significant differences in the expression of glycosylation-related genes in both glioma and non-tumor brain tissues using the heat map and violin plot.Except for LMAN1L, FUT2, SPTB, ADAMTS18, SPTBN5, GALNT12, ST8SIA3, and SPTBN4, all glycosylation-related genes had considerably higher expression levels in tumors than in control tissues (p < 0.001, Fig. 2A,B).

Consensus clustering of glycosylation-related regulators
We used the "ConsensusClusterPlus" package to perform a clustering analysis of the expression matrix of glycosylated genes in each sample to determine the ability of glycosylated genes to discriminate between glioma patients, and the results showed that the area under the Cumulative Distribution Function (CDF) curve with the area under the cluster count k = 2 was considered the best cluster (Fig. 3A-C).We separated the sample into two groups (cluster 1 and cluster 2) and used a PCA 3D pattern plot to show a clear separation between cluster 1 and cluster 2 (Fig. 3D).We ran a survival analysis of glioma patients in clusters 1 and 2 to see if glycosylation genes could differentiate between them.The survival analysis revealed that cluster2 glioma patients had a higher survival rate (p = 1.468e−07) than cluster 1 glioma patients (Fig. 3E), and the clinicopathological characteristics of the patients were also merged with the gene expression heat map (Fig. 3F), patients in clusters 1 and 2 had significant differences in WHO Grade, Histologic, 1p/19q status, MGMTp methylation, and IDH status (p < 0.001).

GO and KEGG enrichment analysis
To investigate the biological differences between the two clusters, we did a difference analysis on the two clusters and screened out the difference genes that met the criteria of adj-P < 0.05 and |Log 2 FC > 1|, then analyzed the difference genes using GO and KEGG.The differential genes were mostly enriched in the "Proteoglycans in cancer," "ECM-receptor interaction," "Focal adhesion," "Cell adhesion molecules," and "Staphylococcus aureus infection" pathways, according to the KEGG analysis.The results of GO analysis showed that differential genes were significantly enriched in the biological processes such as "gliogenesis," "astrocyte differentiation," "collagen fibril organization," "synapse organization" and "positive regulation of cell adhesion".The differential genes were linked to cellular structures like "collagen-containing extracellular matrix," "collagen trimer," and "complex of collagen trimers," as well as molecular functions like "extracellular matrix structural constituent," "glycosaminoglycan binding," and "collagen binding" (Fig. 4A,B, Supplementary Table 3).In addition, we mapped the gene networks that are engaged in each pathway (Fig. 4C-F).

Construction of a glycosylation-related gene prognostic risk model
We used the univariate Cox regression analysis to find 44 genes linked with OS in glioma patients (p < 0.05) to study the influence of glycosylation-related genes on prognosis (Supplementary Fig. 1).We attempted three machine learning methods to screen for differentially expressed genes for constructing prognostic models, and we were unable to select genes common to all three methods for prognostic modeling due to the small number of significantly different prognostic genes (Supplementary Table 4).Finally, Lasso regression analysis and stepwise Cox regression analysis were used to select seven genes linked with OS-related glycosylation for model development, including BGN, C1GALT1C1L, GALNT13, SDC1, SERPINA1, SPTBN5, and TUBA1C, all of which were risk factors (HR > 1) except GALNT13.The sample risk score was obtained based on the formula: risk score = regression coefficient * expression, i.

Evaluation and external validation of risk prognostic models
Patients were separated into two groups for prognostic analysis: low-risk (n = 312) and high-risk (n = 312) based on their median risk score.The survival time of glioma patients declined with rising risk scores, as evidenced in the scatter plot of survival status and the distribution of risk scores, and the survival time of high-risk patients was lower than that of low-risk patients (Fig. 5D).The area under the curve for 1-year, 3-year, and 5-year overall survival was 0.892, 0.909, and 0.876, respectively, indicating that the prediction model had strong predictive power for survival outcomes (Fig. 5E).The Kaplan-Meier curve of the risk-prognosis model revealed that patients in the high-risk group had a considerably lower survival rate (p < 0.001) than those in the low-risk group (p < 0.001, Fig. 5F).We used the CGGA database of glioma patients' gene expression profiles and survival information for   www.nature.com/scientificreports/area under the ROC curve of 0.776 at 1 year, 0.829 at 2 years, and 0.826 at 3 years (Fig. 5G).The Kaplan-Meier survival analysis revealed that high-risk patients had a considerably shorter overall survival (p < 0.001) than low-risk patients (Fig. 5H).The results of combining the CGGA and TCGA-GTEx cohorts revealed that the risk-prognosis model based on glycosylation-related regulators had good predictive potential for glioma patient survival.The expression heat map of the seven glycosylation-associated genes for which the model was built in combination with clinicopathological variables was plotted, and the high-and low-risk groups showed a substantial difference in gene expression (Supplementary Fig. 2).

Model regression analysis and establishment of prognostic nomogram
We used univariate and multivariate Cox analysis to find that age, tumor grade, IDH status, and risk score are all independent prognostic variables for glioma patients (p < 0.05, Fig. 6A,B).The risk score was more accurate than other clinicopathological criteria for prediction, according to multivariate ROC curves (Fig. 6C).We constructed column line graphs of patient mortality at 1, 3, and 5 years to make it easier to forecast glioma patient survival.The nomogram included independent influencing factors such as age, tumor grade, IDH status, and risk score, and the concordance index (C-index) = 0.871 of the column line graphs revealed that risk score had a greater impact on prognosis in glioma patients than WHO Glioma Grade and IDH status (Fig. 6D).The 1-, 3-, and 5-year calibration curves in nomogram revealed good agreement between the nomogram risk model's predicted survival and actual survival (Fig. 6E).

Gene set enrichment analysis
In order to better understand the function of these seven genetic risk profiles, we divided 624 glioma patients into a high-risk group (n = 312) and a low-risk group (n = 312) based on their risk scores and performed GSEA on all genes in the gene expression profiles of the patients in the two groups, respectively (Supplementary Table 5).
The results of the GSEA showed that several immune and tumor-related pathways were involved in the high-risk group, such as JAK-STAT signaling pathway, Antigen processing and presentation, Leukocyte transendothelial migration, B cell activation involved in immune response, integrin mediated signaling pathway and p53 signaling pathway.In addition, these genes are involved in pathways such as glycan biosynthesis and degradation and aminoglycan metabolism, such as N-Glycan Biosynthesis, glycosaminoglycan degradation, amino sugar and nucleotide sugar metabolism, Galactose metabolism and Aminoacyl-tRNA biosynthesis.The low-risk group was mainly enriched in neurotransmitter transmission-related pathways, such as glutamate receptor signaling pathway, GABA receptor binding and Neurotransmitter receptor complexes (Fig. 7A,B).

Immuno-correlation analysis
We performed immune cell and immune function analyses on the risk groups because the risk scores of the seven genes revealed in the enrichment analysis were linked to immune-related pathways that may play an important role in immune regulation.The findings revealed that all immune cells, except eosinophils, had higher immunological scores in the high-risk group of glioma patients than in the low-risk group, including macrophages, natural killer cells, and activated dendritic cells, and so on (Fig. 7C).The immune function evaluation also revealed that glioma patients in the high-risk group had significantly higher immunological scores than those in the low-risk group in a number of immune function modulations, such as Check point, Type I/II IFN Response, and Inflammation promoting (Fig. 7D).Thus, we compared the expression of immune checkpoints in two groups of glioma patients, which are currently undergoing clinical trials in several clinical trials, and the analysis showed that a variety of immune checkpoints were significantly more highly expressed in the high-risk group compared to the low-risk group (p < 0.05), including CD27, CD96, CD274 (PD-L1), CD276 (B7-H3), CTLA4, HAVCR2 (TIM-3), ICOS, IL4I1, LAG3, LGALS9, PDCD1 (PD-1), TIGIT, TNFRSF4 (OX40), TNFRSF9 and TNFRSF18 (GITR) (Fig. 8).

Discussion
More alternatives for treating glioma are being developed all the time, and a combination of microsurgical  prognoses and medication sensitivity, including isocitrate dehydrogenase IDH mutations, 1p/19q co-deletions, and MGMT promoter methylation 16 .IDH mutations and chromosome 1p/19q co-deletions are two of the most prevalent favorable genetic alterations in low-grade gliomas, and patients with IDH mutations and chromosome 1p/19q co-deletions have a better prognosis 14,17 .The addition of the adjuvant drug-alkylating agent temozolomide, which crosses the natural blood-brain barrier to reach intracranial lesions, leading to cellular DNA damage and eventually apoptosis, is the usual treatment for the most malignant glioblastomas; unfortunately, most patients do not benefit from drug therapy due to the presence of the MGMT promoter DNA repair gene 18 .
Existing medicines only improve the prognosis of glioma patients to a limited extent, and additional effective biomarkers for individualized and accurate therapy are urgently needed.
Stable cell surface receptor signaling, cell-matrix interactions, antigen-antibody interactions, cellular infiltration, tumor invasion, and cell motility are all affected by glycosylation modifications 19 .A growing body of evidence suggests that abnormal protein glycosylation is a hallmark of cancer, and that tumor cell surface glycosylation characteristics are one of the key epigenetic changes that occur during the progression of malignant illness 20 .The major glycosylation alterations in gliomas include aberrant N-linked glycosylation and O-linked glycosylation on integrins and receptor tyrosine kinase, as well as aberrant glycoprotein sialylation.Overexpression of glycosyltransferases can promote glycan formation and increase glioma invasion and metastasis, and these glycosyltransferase target genes have great diagnostic and prognostic potential in gliomas 21,22 .We first used consistent clustering to divide glioma patients into two subgroups with optimal k = 2. PCA analysis revealed a clear separation between the two subgroups, so we used GO and KEGG enrichment analysis of differentially expressed genes between the two subgroups, where the KEGG pathway is involved in Proteoglycans in cancer, ECM-receptor interactionFocal adhesion, Cell adhesion molecules, and so on.Positive control of gliogenesis, collagen synthesis, glycans, and cell adhesion related biological processes such as "gliogenesis", "collagen fibril organization", "glycosaminoglycan binding" and "positive regulation of cell adhesion" are all engaged in GO enrichment analysis.These findings point to a significant biological function for differential glycosylation-related genes in the progression of glioma illness.In a previously constructed glioma prognostic model, apoptosis-related indicators showed good predictive value 23 .Currently, studies on the prognosis of gliomas have also shown that annexin A2 (ANXA2) is an unfavorable factor in the prognosis of gliomas, whereas annexin A2 is closely associated with glycosylation, so it is reasonable to believe that glycosylation and related genes play an important role in glioma development 24 .Another paper, using new machine learning algorithms, provides a good evaluation of the predictive efficacy of DNA methylation in the prognosis of glioma patients 25 .This is the first study to use risk profiles of glycosylation-related regulators to predict the prognosis of glioma patients, according to a review of the literature.www.nature.com/scientificreports/44 glycosylation regulators associated with overall glioma survival were identified, and risk prognostic models based on seven genes (BGN, C1GALT1C1L, GALNT13, SDC1, SERPINA1, SPTBN5, TUBA1C) were developed and validated using CGGA database data for external data validation (Fig. 11).In both the prediction model cohort (TCGA-GTEx) and the external validation cohort (CGGA), Kaplan-Meier survival analysis revealed that patient survival was considerably worse in the high-risk group than in the low-risk group.The accuracy of the model was assessed using subject working characteristic curves, with AUCs of 0.892, 0.909, and 0.876 for the TCGA-GTEx cohort predicting 1-, 3-, and 5-year overall survival in glioma patients, respectively, and AUCs of 0.776, 0.829, and 0.826 for the CGGA validation cohort predicting 1-, 3-, and 5-year overall survival in glioma patients.After that, we ran univariate and multifactorial Cox regression analysis, finding that risk score was an independent risk factor for glioma prognosis.To help clinicians visually assess the survival of glioma patients, we used independent prognostic criteria to build a column line graph that could estimate a patient's death at 1, 3, and 5 years, and the calibration curve and C-index revealed that the column line graph was accurate.So far, we've determined that glycosylation-related genes in the model have a good predictive effect on glioma prognosis.We discovered that these seven glycosylation genes are abnormally expressed in a variety of malignancies and are engaged in glycosylation modification-mediated carcinogenesis and immune modulation through a study of the literature.Biglycan (BGN) belongs to a family of tiny leucine-rich proteoglycoproteins that are extracellular matrix components and have a role in skeletal muscle growth and development, collagen fiber assembly, inflammatory modulation, and innate immunity.BGN was overexpressed in pancreatic cancer and inhibited cell development by blocking mitotic G1 phase of pancreatic cancer cells, while BGN was overexpressed in gastric cancer and positively linked with tumor cell repair, invasion, and migration 26,27 .Overexpression of C1GALT1C1 in human colon cancer cells significantly enhances cell migration and invasion, with activation of the EMT signature as the underlying mechanism.C1GALT1C1L is a paralog of C1GALT1C1, a molecular chaperone located in the endoplasmic reticulum that plays a crucial role in proteinO-linked glycosylation 28 .GALNT13 is a member of the N-acetylgalactosaminyltransferase family that is expressed particularly in neuronal cells and is responsible for the synthesis of O-glycan (Tn antigen) 29 .According to studies, GALNT13 is highly expressed in lung cancer and is linked to a bad patient prognosis 30 .SDC1 (syndecans), a member of the type I transmembrane protein family, has higher expression in gliomas than normal brain tissue and is related with tumor aggressiveness and poor prognosis.SDC1 is also related to integrins and plays a role in cell adhesion and migration 31,32 .SERPINA1 belongs to the serine protease inhibitor superfamily, and it has been studied as a target of abnormal www.nature.com/scientificreports/protein fucosylation in pancreatic cancer.Elevated SERPINA1 fucosylation was found to be positively correlated with TNM stage of pancreatic cancer, with high SERPINA1 indicating a poor prognosis 33 .TUBA1C is one of the α-microtubulin isoforms, which has a role in the development of many malignancies.TUBA1C expression was shown to be substantially expressed in low-grade gliomas and to be an independent risk factor for overall survival.TUBA1C expression was also found to be positively connected with the degree of infiltration of various immune cells in low-grade gliomas 34 .In addition, we conducted a GSEA enrichment analysis based on risk scores, with patients in the high-risk group being linked to pathways such as cytokine signaling, inflammatory responses, and immune regulation, as well as glycan synthesis and metabolic function, and those in the low-risk group being linked to synaptic receptors and neurotransmitter transmission.Finally, functional experiments were performed to evaluate the roles of risk genes in the cell viability of glioma cells, which demonstrated that silencing BGN, SDC1, SERPINA1, TUBA1C, C1GALT1C1L and SPTBN5 could inhibit the growth and viability of glioma cells.These findings strengthened the prognostic potentials of our predictive signature in glioma.Malignant cell proliferation and tumor migration are regulated by abnormal glycosylation changes, which also increase tumor-induced immunomodulatory responses 35 .Traditional combination therapy regimens have long been overwhelmed by central system malignancies due to natural anatomical barriers and tumor tissue heterogeneity, and thanks to high-throughput sequencing and the discovery of perivascular glial-lymphatic structures in the central system, immunotherapy and targeted therapy are once again the ultimate weapons clinicians have been waiting for 36 .Immunological checkpoint molecules, tumor-associated macrophages, and dendritic cell vaccines are currently the most appealing biomarkers, and clinical studies of targeted immune combination therapy for glioma are moving in that direction.Treatments targeting the immune checkpoints IDO, CTLA-4, and PD-L1 reduced the number of tumor-infiltrating Tregs cells and improved mouse survival in a glioblast mouse model in previous studies, and combination targeted immunosuppressive therapy has high potential clinical value in high-grade malignant glioma 37 .Microglia, also known as CNS tumor-associated macrophages, are a key component of the glioma tumor microenvironment, and they can produce a range of growth factors and cytokines to control tumor proliferation and cancer cell migration in the glioma immune milieu 38 .Cytokines are important coordinators of the microenvironment in which gliomas develop, and thus cytokine regulators such as Suppressor of Cytokine Signaling 3 (SOCS3) play an active role in neural tissue development and cell proliferation.Thus, glycosylation may influence immune escape from gliomas by regulating the expression levels of certain cytokines in the body 39 .Immune responses against primary and metastatic cancers can be induced by dendritic cell tumor vaccines, and clinical trials have showed remarkable immunotherapeutic potential 6 .As a result, we compared the immune cell infiltration scores and immune function evaluations between the high-risk and low-risk groups, and the results revealed that a variety of immune cell infiltrations, such as Macrophage, activated dendritic cell, Regulatory T cell (Treg), and Natural killer cell, were increased in the high-risk group.In terms of cytolytic activity, checkpoint activation, and inflammation promotion, the immune function of the high-risk group outperformed that of the low-risk group.The expression of the more well-known immune-related targets in glioma, such as CTLA4, PDCD1 (PD-1), CD274 (PD-L1), LAG-3, and CD47, was then examined.Anti-CTLA-4, anti-PD-1, and anti-PD-L1 therapy were found to improve the antitumor immune response of glioblastoma patients treated with targeted CD73.It is suggested that tyrosine metabolizing enzymes regulate adaptive immune processes in gliomas, which are mechanistically related to the remodeling of multiple metabolizing enzymes that induce the expression of programmed death ligand 1 40 .LAG3 is found on activated immune cells and, like PD-1, stimulates tumor cell immune escape.CD47, also known as "integrin-related protein," is thought to enhance the invasion and progression of high-grade glioma.Anti-LAG-3 and anti-CD47 inhibitors could be used as immunotherapy modulators 41,42 .Researchers have proposed multi-immune checkpoint combination therapy and immunotherapy in combination with other targeted pathways in response to the tendency of single-target therapy to relapse and develop drug resistance 43 .These findings led us to believe that risk profiles based on glycosylation modification modifiers linked to glioma tumor formation and immune modulation have a fair chance of accurately predicting glioma patient prognosis.

Figure 1 .
Figure 1.Differential analysis and correlation analysis of glycosylation genes.(A) Volcano plot of differential expression of glycosylation-related genes.Red: up-regulated genes; blue: down-regulated genes; gray: genes with no significant difference in expression.Conditions were met: (FDR) < 0.05 and |Log 2FC|> 1. (B) Pie charts of the proportion of up-and down-regulated genes.Red: up-regulated genes; blue: down-regulated genes.(C) Heat map of correlation between significantly different genes.Red: positive correlation; blue: negative correlation.Numbers in the circles indicate the correlation coefficients between genes.

Figure 2 .Figure 3 .
Figure 2. Differential expression of glycosylation genes in gliomas.(A) Heat map of differential expression of glycosylation-related genes.N: normal tissue; T: glioma tissue; red: high expression; blue: low expression.(B) Differential expression of glycosylation genes in glioma tissues and control tissues violin plot.Red: tumor tissues; blue: normal tissues.

Figure 4 .Figure 5 .
Figure 4. GO and KEGG enrichment analysis.(A) Bar graph of KEGG and GO enrichment analysis with − Log10 (p-adjust) as the horizontal coordinate.(B) Bubble plot of KEGG and GO enrichment analysis with the horizontal coordinate of the enrichment ratio GeneRatio.(C-F) Pathway-related gene network map of Biological Process, Cellular Component, Molecular Function and KEGG Pathway.BP: biological process; CC: cell component; MF: molecular function; GO: gene ontology; KEGG: Kyoto Encyclopedia of Genes and Genomes.

Figure 6 .Figure 7 .
Figure 6.Single-and multi-factor Cox regression analysis and nomogram creation and calibration curves.(A) Forest plot of single-factor Cox regression analysis.(B) Forest plot of multi-factor Cox regression analysis.(C) Multivariate ROC curves.(D) The nomogram incorporating independent prognostic factors age, risk score, WHO Grade and IDH status to predict mortality in glioma patients.(E) The nomogram of 1-year, 3-year and 5-year calibration curves.*p < 0.05, ***p < 0.001.AUC: area under ROC curves.

Figure 11 .
Figure 11.The flow chart demonstrates the main steps of this study.